Entropy of entanglement and correlations induced by a quench: 
Dynamics of a quantum phase transition in the quantum Ising model 
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Quantum Ising model in one dimension is an exactly solvable example of a quantum phase transi- 
tion. We investigate its behavior during a quench caused by a gradual turning off of the transverse 
bias field. The system is then driven at a fixed rate characterized by the quench time tq across the 
critical point from a paramagnetic to ferromagnetic phase. In agreement with Kibble-Zurek mech- 
anism (which recognizes that evolution is approximately adiabatic far away, but becomes approxi- 
mately impulse sufficiently near the critical point), quantum state of the system after the transition 
exhibits a characteristic correlation length £ proportional to the square root of the quench time tq: 
£ = -Jtq. The inverse of this correlation length is known to determine average density of defects 
(e.g. kinks) after the transition. In this paper, we show that this same £ controls the entropy of 
entanglement, e.g. entropy of a block of L spins that are entangled with the rest of the system after 
the transition from the paramagnetic ground state induced by the quench. For large L, this entropy 
saturates at g log 2 £, as might have been expected from the Kibble-Zurek mechanism. Close to the 
critical point, the entropy saturates when the block size L w £, but - in the subsequent evolution 
in the ferromagnetic phase - a somewhat larger length scale I = ^Jtq In tq develops as a result of a 
dephasing process that can be regarded as a quantum analogue of phase ordering, and the entropy 
saturates when L w I. We also study the spin-spin correlation using both analytic methods and real 
time simulations with the Vidal algorithm. We find that at an instant when quench is crossing the 
critical point, ferromagnetic correlations decay exponentially with the dynamical correlation length 
£, but (as for entropy of entanglement) in the following evolution length scale I gradually develops. 
The correlation function becomes oscillatory at distances less than this scale. However, both the 
wavelength and the correlation length of these oscillations are still determined by £. We also derive 
probability distribution for the number of kinks in a finite spin chain after the transition. 

PACS numbers: 03.65.-w, 73.43.Nq, 03.75.Lm, 32.80.Bx, 05.70.Fh 



INTRODUCTION 



Phase transition is a fundamental change in the char- 
acter of the state of a system when one of its parameters 
passes through the critical point. States on the opposite 
sides of the critical point are characterized by different 
types of ordering. In a second order phase transition, the 
fundamental change is continuous and the critical point 
is characterized by divergences in the coherence (or heal- 
ing) length and in the relaxation time. This critical slow- 
ing down implies that no matter how slowly a system is 
driven through the transition its evolution cannot be adi- 
abatic close to the critical point. If it were adiabatic, then 
the system would continuously evolve between the two 
types of ordering. However, in the wake of the necessar- 
ily non-adiabatic (and approximately impulse) evolution 
in the critical region, ordering of the state after the tran- 
sition is not perfect: It is a mosaic of ordered domains 
whose finite size depends on the rate of the transition. 
This scenario was first described in the cosmological set- 
ting by Kibble |l[ who appealed to relativistic casuality 
to set the size of the domains. The dynamical mecha- 
nism relevant for second order phase transitions was pro- 
posed by one of us 2]. It is based on the universality 
of critical slowing down, and leads to prediction that the 
size of the ordered domains scales with the transition 



time tq as Tq , where w is a combination of critical expo- 
nents. The Kibble-Zurek mechanism (KZM) for second 
order thermodynamic phase transitions was confirmed by 
numerical simulations of the time-dependent Ginzburg- 
Landau model @ and successfully tested by experiments 
in liquid crystals [4k superfluid helium 3 [5], both high- 
T c Q and low-T c ^superconductors and even in non- 
equilibrium systems [8] . With the exception of superfluid 
4 He - where the early detection of copious defect forma- 
tion [9( was subsequently attributed to vorticity inadver- 
tently introduced by stirring (To| . and the situation re- 
mains unclear - experimental results are consistent with 
KZM, although more experimental work is clearly needed 
to allow for more stringent experimental tests of KZM. 

The Kibble-Zurek mechanism is thus a universal the- 
ory of the dynamics of second order phase transitions 
whose applications range from the low temperature Bose- 
Einstein condensation (BEC) [ll| to the ultra high tem- 
perature transitions in the grand unified theories of high 
energy physics. However, the zero temperature quan- 
tum limit remained unexplored until very recently, see 
Refs.JH, 0, Il4l [TEl and an example of disordered 
system in Ref. |l7| . and quantum phase transitions are 
in many respects qualitatively different from transitions 
at finite temperature. Most importantly time evolution is 
unitary, so there is no damping, and there are no thermal 
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fluctuations that initiate symmetry breaking in KZM. 

Quantum state of the many body system is indeed pro- 
foundly different from a classical state: Instead of a single 
broken symmetry configuration it may (and, generally, 
will) contain all of the possible configurations in a super- 
position. In addition to the 'classical' way of character- 
izing this state through the density of excitations (e.g. 
defects), one can wonder how entangled various parts of 
the system are with each other. Von Neumann entropy 
of a fragment due to its entanglement with the rest of 
the system is a convenient way to quantify this. It can 
be computed as a function of the size of the fragment. 
In equilibrium, and away from the critical point, this en- 
tropy of entanglement saturates at distances of the order 
of the coherence length £ of the system at values ~ In £ 
for one dimensional systems [181 ]. However, at the criti- 
cal point (where equilibrium coherence length £ becomes 
infinite) entropy of entanglement diverges with the size 
of the fragment. In particular, in one dimensional sys- 
tems, the critical entanglement entropy diverges logarith- 
mically (~ In L) with the length L of the chain fragment 

Unlaid. 

This equilibrium behavior suggests the question: What 
is the entanglement entropy left behind by an out-of- 
equilibrium phase transition? Such a transition will pass 
through the critical point (where entanglement entropy is 
logarithmically divergent) but this will happen at a finite 
rate set by the quench time tq . We show that the result- 
ing entanglement entropy is of the order of ~ ln£, where 
£ is the healing (coherence) length at the instant when 
critical slowing down forces the system to switch from 
the approximately adiabatic to approximately impulse 
('diabatic') behavior. This suggests that the same pro- 
cess that determines the size of regions that "break sym- 
metry in unison" (which sets the density of topological 
defects left by the transition) is also responsible for the 
resulting entanglement of formation left by the quench. 
This finding is consistent with recent results on quan- 
tum phase transitions induced by instantaneous quenches 
[2H I22I [23j which indicate that structures present in the 
initial pre- transit ion state determine the structures (and 
hence entanglement of formation) that arise after an in- 
stantaneous quench: Our results also suggest that - in 
accord with KZM - it is a good approximation to con- 
sider quench to be approximately adiabatic until the in- 
stant t r~j 1/-Jtq before the critical point is reached, and 
approximately impulse (e.g. nearly instantaneous) inside 
this time interval. This also confirms and extends re- 
sults of the recent study of Cherng and Levitov [l6[ who 
computed entropy density and correlations induced by 
quenches in one-dimensional chains, and concluded that 
their results support KZM. 

While our results below are established for the one- 
dimensional quantum Ising model (which has the consid- 
erable advantage of being exactly solvable), we conjecture 
that similar behavior will be encountered in other quan- 
tum phase transitions, and that their non-equilibrium 
evolution can be anticipated using equilibrium critical ex- 



ponents using KZM. This conjecture can be then tested 
in a variety of systems that undergo quantum phase tran- 
sitions both in condensed matter and in atomic physics 
experiments. 

According to Sachdev [24| , the understanding of quan- 
tum phase transitions is based on two prototypical mod- 
els. One is the quantum rotor model and the other is the 
one-dimensional quantum Ising model. Of the two only 
the Ising model is exactly solvable. It is defined by the 
Hamiltonian 

JV 

H = -EG?<+«+i) • (!) 

71=1 

with periodic boundary conditions 

= <?i . (2) 

Quantum phase transition takes place at the critical value 
g = 1 of an external magnetic field. When g ^> 1, the 
ground state is a paramagnet | -^^^ ■ ■ ■ — >) with all 
spins polarized along the x-axis. On the other hand, 
when g -C 1, then there are two degenerate ferromagnetic 
ground states with all spins pointing either up or down 
along the z-axis: | TTT - - - T> or I III • ■ ■ I)- In an h> 
finitesimally slow classical transition from paramagnet to 
ferromagnet, the system would choose one of the two fer- 
romagnetic states. In the analogous quantum case, any 
superposition of these two states is also a 'legal' ground 
state providing it is consistent with other quantum num- 
bers (e.g. parity) conserved by the transition from the 
initial paramagnetic state. However, when N — ► 00, then 
energy gap at g = 1 tends to zero (quantum version of 
the critical slowing down) and it is impossible to pass 
the critical point at a finite speed without exciting the 
system. As a result, the system ends in a quantum su- 
perposition of states like 

I ... TIUUTTTTTTTIUITTTTTTI •••) (3) 

with finite domains of spins pointing up or down and sep- 
arated by kinks where the polarization of spins changes 
its orientation. Average size of the domains or, equiva- 
lently, average density of kinks depends on a transition 
rate. When the transition is slow, then the domain size 
is large, but when it is very fast, then orientation of indi- 
vidual spins can become random, uncorrelated with their 
nearest neighbors. Transition time tq can be unambigu- 
ously defined when we assume that close to the critical 
point at g = 1 time-dependent field g(t) driving the tran- 
sition can be approximated by a linear quench 

g(t < 0) = — . (4) 

with the adjustable quench time tq. Density of kinks af- 
ter the linear quench was estimated in Ref. [l3| showing 
that KZM can be also applied to quantum phase transi- 
tions. In this derivation, it is convenient to use instead of 
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](t) a dimensionless parameter e(t) 



1. As 



in classical transitions |2j, one can assume the adiabatic- 
impulse approximation [2^, |26| . The quench begins in 
the ground state at large initial g and the initial part 
of the evolution is adiabatic: the state follows the in- 
stantaneous ground state of the system. The evolution 
becomes non-adiabatic close to the critical point when 
the energy gap — |e| becomes comparable to the instan- 
taneous transition rate |e/e|. This condition leads to an 
equation solved by t = ^Jtq - the instant when the adia- 
batic to impulse transition occurs - which in turn yields 

A —1/2 

e ~ Tq and corresponds to the coherence length in the 



ground state @, [HI 



r l/2 



(•5) 



Assuming impulse approximation, the quantum state 
does not change during the following non-adiabatic stage 
of the evolution between e and — e. Consequently, the 
quantum state at — e is expected to be approximately the 
ground state at e with the coherence length proportional 
to £ and this is the initial state for the final adiabatic 
stage of the evolution after e. This argument shows that 
when passing across the critical point, the state of the 
system gets imprinted with a finite KZ correlation length 
proportional to £. In particular, this coherence length de- 
termines average density of kinks after the transition as 



1 



.1/2 



(6) 



Q 



This is an order of magnitude estimate with an unknown 
0(1) prefactor. The estimate was verified by numerical 
simulations in Ref. • Not much later the problem was 
solved exactly in Ref. see also Ref. jig , with the 
exact solution confirming the KZM scaling in Eq. ([6]). 

In the next section we review and expand the exact 
solution, and then use the expanded version to obtain a 
more complete set of results. In Subsection III D[ we de- 
rive Gaussian probability distribution for the number of 
kinks measured after a quench in a finite Ising spin chain. 
In Section UTTl we calculate entropy of entanglement of a 
block of L spins after a dynamical transition and in Sec- 
tion [IV] we work out spin-spin correlation functions. We 
conclude in Section IVl 



II. EXACT SOLUTION 
A. Energy spectrum 

Here we assume that the number of spins N is even for 
convenience. After the nonlocal Jordan- Wigner transfor- 
mation m, 



at 



1 



+ c n) n( i - 2c ™ c ™) . 



(7) 

(8) 



introducing fermionic operators c n which satisfy anti- 
commutation relations {c m ,4} = S mn and {c m ,c„} = 
{cj^c^} =0 the Hamiltonian (fT]) becomes [H[ 



H 



p + H + p + + p - H - p - 



(9) 



Above; 



JY 



n=l 



N 



1 ± JJ (1 - 24c) 



71=1 



(10) 

are projectors on the subspaces with even (+) and odd 
(— ) numbers of c-quasiparticles and 



iV 



9 c n c n c n c n -|_ i c n i c n 



§ + >,. 



(11) 



are the corresponding reduced Hamiltonians. The c„'s in 
H~ satisfy periodic boundary conditions cjv+i = c\, but 
the Cji's in H + must obey cn+i = — Ci - e.g, "antiperi- 
odic" boundary conditions. 

The parity of the number of c-quasiparticles is a good 
quantum number and the ground state has even parity 
for any non-zero value of g. Assuming that a quench be- 
gins in the ground state, we can confine to the subspace of 
even parity. H + is diagonalized by the Fourier transform 
followed by Bogoliubov transformation [28| . The Fourier 
transform consistent with the antiperiodic boundary con- 
dition cj\r + i = — ci is 



-17T/4 



E 



c k e 



(12) 



where pseudomomenta k take "half- integer" values: 



± 



1 2vr 



N - 1 2vr 



2 N' 2 
It transforms the Hamiltonian into 

H+ = { 2 l5 " cos(fca)]4c fe - 



N 



(13) 



sin(fca) 



C k C -k + c -kCk 
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(14) 



Diagonalization of H + is completed by the Bogoliubov 
transformation 



cfe = u k jk + v*_ k ^_ k , 



(15) 



provided that Bogoliubov modes (uk,Vk) are eigenstates 
of the stationary Bogoliubov-de Gennes equations 

e Uk = +2[<? — cosk]uk + 2sinkvk , 

e Vk — — 2[g — cos k]vk + 2 sin kuk ■ (16) 

There are two eigenstates for each k with eigenenergies 
e = ±6^, where 



eft = 2 \J [g — cos k] 2 + sin 2 k . 



(17) 
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The positive energy eigenstate 



(u k , Vk) ~ (g — cos fc) + \/ g 1 — 2g cos k + 1, sin fc 



which has to be normalized so that |iifc| 2 + \vk\ 



(18) 
= 1, 



defines the quasiparticle operator jk = u k c k + w -fcCl fc , 
and the negative energy eigenstate (u k , vj~) = (—Vk, u k ) 



defines 



7* 



Cfc 



-k c -k 



After 



the Bogoliubov transformation, the Hamiltonian i7 + 



I Efc £ fc (7fc' 



7fc7fe - Ik ik ) is equivalent to 



Eft 



7 fc 7fc 



(19) 



This is a simple-looking sum of quasiparticles with half- 
integer pseudomomenta. However, thanks to the projec- 
tion P + H + P + in Eq. © only states with even numbers 
of quasiparticles belong to the spectrum of H. 



B. Linear quench 



ihf t c u 



[ck,H + ] with the constraint that ■^•jk = 0. 
The Heisenberg equation is equivalent to the dynamical 
version of the Bogoliubov-de Gennes equations (fT6|) : 



ih—Uk = +2 [g(t) — cosfc] Uk + 2 sink Vk , 

ifi.^Wfc = — 2 [g(t) — cosfc] Vk + 2sinfc Uk ■ (23) 

At any value of g, Eqs. (|23j) have two instantaneous eigen- 
states. Initially, the mode [iik{t), Vk{t)] is the positive en- 
ergy eigenstate, but during the quench it gets "excited" 
to a combination of the positive and negative mode. At 
the end of the quench at t — when g = we have 



[u k (0),v k (0)} = A k (u k ,v k ) 



B k {u- k ^) (24) 



and consequently j k = A k j k — B k j k . The final state 
which is, by construction, annihilated by both "f k and 

7-fc is 



I^(0)) = n(^ + ^7fc7ifc)|0) 



(25) 



fc>0 



In the linear quench Eq. ([5]), the system is initially 
(t <C — tq) in its ground state at large initial value of 
g > 1, but when g is ramped down to zero, the system 
gets excited from its instantaneous ground state and, in 
general, its final state at t = has finite number of kinks. 
Comparing the Ising Hamiltonian Eq. |T]) at g — with 
the Bogoliubov Hamiltonian (JT9J) at g = we obtain a 
simple expression for the operator of the number of kinks 



1 N 

-Y 

2 ^ 



(1 



H 7 fc7fc 
fc 



(20) 



The number of kinks is equal to the number of quasipar- 
ticles excited at g = 0. The excitation probability 



Pk - (m\ihk\m) 



(21) 



in the final state can be found with the time-dependent 
Bogoliubov method. 

The initial ground state is Bogoliubov vacuum |0) an- 
nihilated by all quasiparticle operators "f k which are de- 
termined by the asymptotic form of the (positive en- 
ergy) Bogoliubov modes (u kl v k ) fa (1,0) in the regime 
of g 1. When g(t) is ramped down, the quantum 
state \ip{t)} gets excited from the instantaneous ground 
state. The time-dependent Bogoliubov method makes an 
Ansatz that \tp(t)) is a Bogoliubov vacuum annihilated 
by a set of quasiparticle annihilation operators "/ k defined 
by a time-dependent Bogoliubov transformation 



Cfc = u k (t)% + v*_ k (t)jl_ 



(22) 



with the initial condition [u k (— oo), v k (— oo)] = (1,0). 
In the Heisenberg picture, the Bogoliubov modes 
[u k (t),v k (t)] must satisfy Heisenberg equation 



Pairs of quasiparticles with pseudomomenta (fc, 
excited with probability 



Pk 



\B k \ 2 



-fc) are 



(26) 



which can be found by mapping Eqs. (|23[) to the Landau- 
Zener (LZ) problem (similarity between KZM and LZ 
problem was first pointed out by Damski in Rcf. [25]). 
The transformation 



r.^-nM-Uo.* 



brings Eqs. (|23[) to the standard LZ form [2£ 

d 1 / A \ 1 

ih—Uk = --{rAkjuk + -v k , 

(XT i I 

• t d 1, A x 1 

in—vk = +-{TA k )Vk + -Uk , 
or i i 



(27) 



(28) 



with A,: 1 



4rn sin fc. Here the time r runs from — oo 



to rg na j = 2tq sin(2fc) corresponding to t — 0. Tunnel- 
ing between the positive and negative energy eigenstates 
happens when r G (— A^AJT 1 ). Tfi na i is well outside 
this interval, Tfi na i ^> A^ 1 , for long wavelength modes 
with \k\ <C ?. For these modes, time r in Eqs. (|28|) can 
be extended to +00 making them fully equivalent to LZ 
equations [29]. This equivalence can be used to easily 
obtain several simple results [15| | described in the next 
subsection. 



C. Simple results 

In the limit of slow transitions we can assume that only 
long wavelength modes, which have small gaps at their 
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anti-crossing points, can get excited. For these modes, we 
can use the LZ formula [29j for excitation probability: 



Pk 



-1-KTQk 1 



(29) 



This approximation is self-consistent only when the 
width of the obtained Gaussian (47ttq) -1 / 2 is much less 
than j or, equivalently, for slow enough quenches with 
tq 3> 1. With the LZ formula (f2"9"| . we can calculate the 
number of kinks in Eq. ([20]) as 



TV = 



k 



Pk 



(30) 



• When N — > oo the sum in Eq. ([30)1 can be replaced 
by an integral. The expectation value of density of 
kinks becomes 



There are at least two interesting special cases: 



r N 1 

n = lim — = — 



dk p k 



1 1 



(31) 



— 1/2 

The density scales like Tq in agreement with 
KZM, see Eq. (|6|). Thus, slower quenches lead to 
fewer defects. 

• Following Ref. [l3| , we can ask what the fastest tq 
is when no kinks get excited in a finite chain of size 
N . This critical tq marks a crossover between adi- 
abatic and non-adiabatic regimes. In other words, 
we can ask what is the probability for a finite chain 
is to stay in the ground state. As different pairs 
of quasiparticles (k,—k) evolve independently, the 
probability to stay in the ground state is the prod- 
uct 



(32) 



fc>0 



Well on the adiabatic side only the pair (jj, — j?) 
is likely to get excited and we can approximate 



!- ex p(- 27r3 ^) ■ (33) 



A quench in a finite chain is adiabatic when 



N 2 



TQ » 2^ 



(34) 



Reading this inequality from right to left, the size 

1 /2 

N of a defect-free chain grows like Tq . This is 
consistent with Eq. (|6l3ip . 



D. Probability distribution of the number of kinks 



Equation ()31|) gives an average density of kinks mea- 
sured after a quench to zero magnetic field g = 0. In 
a finite chain, an average number of kinks is AT = Nn, 



provided that the transition is non-adiabatic unlike in 
Eq. (|34p. The average AT is an expectation value of a 
probability distribution P(Af) for the number of kinks Af 
measured after a quench. 

The number of kinks Af is the number of quasiparti- 
cles excited by the end of the quench. As quasiparticles 
are excited in pairs with opposite quasimomenta [k, —k), 
the number of kinks Af must be even. A pair (k, —k) 
is excited with the probability pk in Eq. (|29|) . We can 
asign to each pair of quasiparticles a random variable 
Xk which is 1 when the pair is excited and otherwise: 
Xk = 1 with probability pk and Xk = with probability 
1 — pk- We want a probability distribution for the sum 
k>0 2xk- This is a sum of independent random 
variables of finite variance so for a large number of vari- 
ables, the sum Af becomes a Gaussian random variable 
with a mean Af — Nn and finite variance ~ N. 

We have to be careful to specify when the number of 
random variables is large. Naively, on a TV-site lattice, 
there are N/2 pairs of quasiparticles [k, —k), or indepen- 
dent random variables x k , and the number seems to be 
large when N is large. On second thought, it is clear that 
variables with pk ~ or pk ~ 1 cannot really count be- 
cause they are hardly random at all. A look at the Gaus- 
sian pk in Eq. (f2"9"| shows that the range of k > where 
<C Pk *C 1 has width ~ —^= which accommodates 

~ — = discrete values of pseudomomentum k. The rele- 

vant number of random variables is ~ -£= ~ Nn = Af. 
It is large when the average number of kinks is large, 



Af > 1 



(35) 



With this assumption P(Af) is Gaussian. 

Keeping this assumption in mind we can proceed as 

P(Af) = £ <W,£ fe>o2 *, 

]J Po,x fc (1 - Pk) + Si,x„Pk] = (36) 

fc>0 

±- f dqe-^J[[{l-p k )+p k e^] = 

71 J- 71 k>0 



j_ r 

2tt 



dqe~ lqN 'P(q) 



(37) 



It is convenient to evaluate first 

\nP(q) = ^ln[(l-p fe )+p fe e 2i «] w 

k>0 

— dk \n[{l-p k )+p k e 2lq ] . (38) 

Here we used N 3> 1, a necessary condition for our as- 
sumption that Af 3> 1. After changing the integration 
variable to u = t^—, we can extend the integration over 

2-KTL 1 ° 
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u to infinity: 

lnP(g) = 77 du In (1 
Jo L 



AT 



V2-1 2 



i 2 ^j _|_ e —lTU 2 +2iq 

0(q 3 ) . 



(39) 



Here we used again AT = Nn. 

Finally, combination of Eqs. (J2U [2Z|) gives the ex- 
pected Gaussian probability distribution 



P(A0 = 
with a variance 



: exp 



M 



a M = 



(2-V2) 77 



(40) 



(41) 



In the course of our approximations, we have lost the 
constraint that M must be even, but this is not a big 
mistake when AT ^> 1. 

It is also interesting to study the opposite adiabatic 
regime when Nn -C 1. When tq is large we can approx- 
imate the product in Eq. ([36]) by a single factor with the 
lowest fc = ir/N, 

P{JV) « ^ <W,2s [&0,a (l - Ftt/Jv) + h,sPir/N\ = 
s=0,l 

<W,0 (l - Piv/n) + <W,2 Ptt/N • (42) 

This is a good approximation when tq 3> 



Eq. (|34|) . In this adiabatic regime 

77 = 2 Pv/N = exp(-27r 3 ^|) 



AT 

27? 



as m 



(43) 



The average number of kinks decays exponentially with 
tq and the KZ power law scaling (|3"Tj) does not extend to 
this adiabatic regime, as was already noted in Ref. 



E. Exact solution and the two scales of length 

So far we have avoided writing down solutions of the 
Landau-Zener equations (l28]l whose general form is, see 
e.g. Appendix B in Ref. [261 ], 



Ufc(r) = - [aD- s -i(-iz) + &L>_ s _i(iz)] , 
9 



Uk{r) 



-A fc r + 2?; 



c9 7 



Vk{r) , 



(44) 



with arbitrary complex parameters a,b. Here D m (x) 
is a Weber function, s = .\ , and z = 



i fc re' 



r/4 



The parameters a, 6 are fixed by the initial conditions 
Uk(— oo) = 1 and Ufc(— oo) = 0. Using the asymptotes of 
the Weber functions when r — > — oo, we get a = and 



-7r/8A fc 

4A fc 



(45) 



The solution of the linear quench problem is then 

Ufc(r) = -bD- s -i(iz) , 

Wfc(r) = (-A k T + 2i-^-j v k (r) , 



(46) 



At the end of the quench for t = and when r = = 
2TQsin(2/c), the argument of the Weber function iz = 
V /A7re i7r / 4 = 2 v /rv2e" r/4 cos(fc)sign(fc). In the limit of 
large tq the modulus of this argument is large for most 
k, except the neighborhoods of fc = ±~, and we can again 
use the asymptotes of the Weber functions. After some 
work we get the products 



\Uk\ 

kl 2 

u k v* k 



1 — cos fc 
2 

1 — |^fc| 2 
— sin k + 



3 — Itttq sin k 



sign(fc) e - nTQSin Vl -e-^Q sin2fc e^ fc 
7T ( A fe r^ : lnA fe | lnr fc 
4 



arg 



2 

rfi 



4A fc 2A fc 



4A* 



(47) 



Here T{x) is the gamma function. 

We expect that for large tq only modes with small 
|fc| <C \ get excited. In this long wave length limit, the 
products can be further simplified to 



kl 2 


1 — cos fc 
= 2 4 


■ e _27rrQ ' c2 


kl 2 


= i-kl 2 , 




u k v* k 


= — sin fc + si 


gn(fc) e^Vl-e- 




7T 

= I +2 ^- 


(2-ln4)T Q fc 2 + 




k 2 r Q In Tq 


-arg[r(l + zr Q fc 2 )] . 



(48) 



These products depend on fc and tq through two combi- 
nations: tq/c 2 , which implies the usual KZM coherence 
length £ = ^/tq, and k 2 TQ In tq which implies a second 

length scale U tq In tq . The final quantum state at g = 
cannot be characterized by a single scale of length. Phys- 
ically, this appears to reflect a combination of two pro- 
cesses: KZM that sets up initial post-transition state of 
the system, and the subsequent evolution that can be 
regarded as quantum phase ordering. 



III. ENTROPY OF A BLOCK OF SPINS 

Von Neumann entropy of a block of L spins due to its 
entanglement with the rest of the system; 



S(L) 



Tr p L log 2 p L , 



(49) 
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is a convenient measure of entanglement. Above is re- 
duced density matrix of the subsystem of L spins. In re- 
cent years this entropy was studied extensively in ground 
states of quantum critical systems [H, E3 EH E2, EE EH, 
I32I . l33l ]. At a quantum critical point, the entropy di- 
verges like log L for large L with a prefactor determined 
by the central charge of the relevant conformal field the- 
ory [lil. Il9l. |20| . In particular, in the quantum Ising model 
at the critical g = 1 



S^(L) ~ -log 2 L 
6 



(50) 



and 



1 



\ Cm Cr>, / — 



dk Ukvl e 



ik(m-n) 



1 r lnr Q >l 

— / dk Ukv'k sinfc(m — n) 



(54) 



sign(m — n) x 



jh,\w.-n\ 



it 2tq- 



[m-nT 



-e 4 ^ 



1 — e 



2\/7T £ Z 



for large L. Slightly away from the critical poin t, the 
entropy saturates at a finite asymptotic value jl8l | 



r-GS 
■^00 



log 2 £ 



(51) 



when the block size L exceeds the finite correlation length 
£ in the ground state of the system. 



A. Entropy after dynamical transition 

In a dynamical quantum phase transition the quantum 
state of the system developes a finite correlation length 
£ ~ ^Jtq. If this dynamical correlation length were the 
only relevant scale of length, then one could expect en- 
tropy of entanglement after a dynamical transition given 
by Eq. (|5ip with £ simply replaced by £. However, as 
we saw in Eq. (|4"5)) . there are two scales of length, and - 
strictly speaking - there is no reason to expect that ei- 
ther of them alone is relevant in general. This is why we 
shall not rely on scaling arguments alone and will go on 
to calculate the entropy of entanglement "from scratch" . 



We proceed in a similar way as in Ref . [20J, |21|, |30j, l32j, 
[33j and define a correlator matrix for the block of L spins 



n = 



, l-a 



(52) 



where a and (3 are Lx L matrices of quadratic correlators 



2tt 



dk \uk\ e 



2 ik(m-n) 



TQ»1 



2°0,|m-n| — J°l,|rre-n| 



T 



e 8 " « 
2\/27r £ 



(53) 



Here £ 
and 



tq is the KZM dynamical correlation length 



tq In TQ 



(55) 



We note that a and (3 are Toeplitz matrices with con- 
stant diagonals. The expectation values (...) are taken 
in the dynamical Bogoliubov vacuum state. As this state 
is Gaussian, all higher order correlators can be expressed 
by the matrices a and - they provide complete char- 
acterization of the quantum state after the dynamical 
transition. The matrices depend on both scales £ and I 
and both scales are necessary to characterize the Gaus- 
sian state. 

As observed in Ref. [13, EH, IM GH El, the entropy 
can be conveniently calculated as 

S(L, tq) = - Tr p log 2 p = - Tr n log 2 n . (56) 



In this calculation we use Eq. (|47J) and Eqs. ([531 \§M 
but without their large tq approximations. The calcu- 
lation involves a numerical evaluation of the integrals in 
Eqs. (|53l54p and numerical diagonalization of the ma- 
trix n. Results are shown in Panel A of Fig[TJ The 
entropy grows with the block size L and saturates at 
a finite value Soo (tq) for large enough L. In Panel B, 
we fit the asymptotic entropy with the linear function 
Soo (tq ) = A + B In tq . The simple replacement of £ by 
£ = ^Jtq in Eq. (|5Tj) suggests the asymptotic value 



0.120 InTQ . (57) 



5*00 (t q ) ~ ^log 2 ^ ~ ^ In r Q 



Our best fit gives B = 0.128 ±0.004 and A = 1.80 ±0.05. 
The best B is in reasonably good agreement with the 
expected value of 0.120. 

In Panels C and D of the same figure, we rescale values 
of entropy S(L,tq) by its asymptotic value Soo(tq) ~ 
A + B\utq. After this transformation we can better 
focus on how the entropy depends on the block size L. 
A simple hypothesis would be that entropy depends on 
£ and saturates when L > £. To check if this is true, in 
Panel C we also rescale the block size L by £ = ^Jtq and 
find that while this rescaling brings plots close to overlap, 
they do not overlap as well as one might have hoped. 
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FIG. 1: Panel A shows entropy of a block of L spins after 
the dynamical phase transition as a function of the block size 
L. The multiple plots correspond to different values of the 
quench time tq . For all tq , the entropy grows with the block 
size L and saturates at a finite value Soo (tq) for large enough 
L. In Panel B, we fit this asymptotic value of entropy with 
the function Soo (tq) — A + Bin tq. The best fit has B = 
0.128 ± 0.004 and A = 1.80 ± 0.05. This B is in reasonably 
good agreement with the expected value of B = ^ = 0.120. 
In Panels C and D, we rescale values of entropy S(L,tq) by 
the best fit to its asymptotic value Soo(tq) = A + Bin tq. 
With this rescaling we can focus on how the entropy depends 
on the block size L. In Panel C, we also rescale the block size 
by C — \/tq an d fi n d that the rescaled plots do not overlap 
exactly. However, as shown in Panel D, rescaling the block 
size L by the second scale I — ^Jtq In tq makes the six plots 
overlap quite well. 



By contrast, as shown in Panel D, rescaling of the block 
size L by I — ^Jtq In tq makes the multiple plots overlap 
quite well indeed. In conclusion, our results support the 
statement that the entropy saturates at 



when 



Soo(T Q ) ~ - log 2 



(58) 



(59) 



i.e. the entropy of a large block of spins is determined 
by Kibble-Zurek dynamical correlation length £ = Jtq, 
but the entropy saturates when the block size is greater 
than the second scale I — ^tqItitq. 

We believe that the KZ correlation length £ is deter- 
mined when the system is crossing the critical point while 
the second longer scale builds up after the system gets 
excited from its adiabatic ground state near the critical 
value of a magnetic field. At the origin of the second scale 
is the non-trivial dispersion relation of excited quasipar- 
ticles. The £;-dependent Ck in Eq. (fTT)) leads to a gradual 
evolution of matrix elements of the correlator II which 



FIG. 2: Panel A shows entropy of a block of L spins during 
the dynamical phase transition at the critical point g = 1 as 
a function of the block size L. The multiple plots correspond 
to different values of the quench time tq. For all the quench 
times, the entropy grows with the block size L and saturates 
at a finite value S x (tq) for large enough L. In Panel B we fit 
this asymptotic value of entropy with the function Soo(tq) = 
A + B In tq. The best fit has B = 0.126 ± 0.005 and A = 
0.80 ±0.03. This B is in reasonably good agreement with the 
expected value of B = = 0.120. In Panels C and D, we 
rescale the entropy S(L, tq) by the best fit to its asymptotic 
value 5*00 (tq) w A + B In tq. With this rescaling we can focus 
on how the entropy depends on the block size L. In Panel D, 
rescaling the block size L by £ = ^/tq makes the six plots 
overlap quite well. 



are given by integrals over k in Eqs. (|53I54|) . To support 
this scenario we calculated the entropy of entanglement 
at the moment when the system is crossing the critical 
point at g = 1 . The results collected in Figure are con- 
sistent with our expectation that near the critical point, 
when the scale I set up by quantum phase ordering only 
begins to build up, £ is still the only relevant scale of 
length. 



B. Impurity of the state after transition 

We were not able to do a fully analytic calculation of 
entropy This is why it may be worthwhile to calculate 
analytically another more easily tractable entanglement- 
related quantity. For example, the "impurity" of the cor- 
relator matrix II 



/(it) = Tr n (1 - n) 



(60) 



is zero only when the L spins are in a pure state i.e. 
when all eigenvalues of II are either or 1 . It is maximal 
when all the eigenvalues are h , or when the state is most 
entangled. Thanks to its simple quadratic form, it can 
be calculated relatively easily. 
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Simple calculation using the block structure of II in 
Eq. (f52"|) and the Toeplitz property of the block matrices 
a and j3 leads to 

j=L-l 

I(L) = 2(La - £ (L-\j\)(a 2 + \P j \ 2 )],{61) 

j=l-L 

where aj — a>j,o and f3j — /3j.Q. aj and f3j can be ex- 
pressed by the inverse Fourier transforms in Eqs. (|53|54| . 
Using normalization \u k \ 2 + \v k \ 2 = 1 and completeness 
of the Fourier basis we notice that 



1 

a Q = — dk \u k \ 2 = 

^ r dk{\ Uk \±+\u k vt\ 2 ) 



(62) 



E (I&i 2 + ki 2 )- 



This leads to 



L-l 



j=L j=l 

Since we assume that tq > 1, we can approximate these 
sums with integrals. Further calculation gives: 



I(L) 



1 1 



Ml!2)fl_ e -i^- 
2tt 2 I 



1 



L 



Erfc [ ^ ) + 
V2Erfc ( ^ ^ - Erfc ( 



(64) 



Here Erfc(cc) is the complementary error function defined 
as: 



Erfc(x) = —= 



' dt. 



(65) 



This impurity saturates at 1^ ~ In rg when the block 
size L ^> I, or in short 



Jo 



2tt 2 



(66) 



The impurity saturates at the second scale I = ^Jtq In tq 
in consistency with our results for the entropy. 

It is interesting to compare the dynamical impurity 
(I66p with the impurity in the ground state of the system. 
Simple calculation gives the asymptote of impurity in the 
ground state at the critical point 



I^(L) 



InL 

tt2 



(67) 



when InL ^ 1. Near the critical point, the asymptote is 
valid for the block size L much less than the correlation 
function L <C £ and at larger L the impurity saturates at 



GS 



(68) 



Simple replacement of £ in this equation by the dynami- 
cal KZ correlation length £ = ^tq gives the asymptotic 
value of the dynamical impurity in Eq. (|66[) . Again, this 
"replacement rule" is the same as for the entropy. 



IV. 



CORRELATION FUNCTIONS 



Correlation functions are of fundamental interest in 
phase transitions because they provide direct manifesta- 
tion of their universal properties and are in general easily 
accesible experimentally. In this Section we present our 
results for spin-spin correlation functions during a dy- 
namical quantum phase transition. 

To begin with, we observe that for symmetry reasons 
the magnetization (er 2 ) = 0, but the transverse magneti- 
zation 



(1 



Hen) 



2oo-l 



1 



2tt V / 2^ 



(69) 



which is valid when tq > 1. This is what remains of the 
initial magnetization (er^) = 1 in the initial ground state 
at g — > oo. As expected, when the linear quench is slow, 
then the final magnetization decays towards (a*) = 
characteristic of the ground state at the final g = 0. 
Final transverse spin-spin correlation function at g = 

is 



4(|/3«| 2 -M 2 ) » 

_7rR 2 / _ tvR 2 

2 i' 2 1 1 — e <i~7^ 



e 2 



7T £ I 



2tt 2 £ 2 



(70) 



(71) 



when R > 1 and In tq 1. This correlation function 
depends on both £ and I. Long range correlations 



(72) 



decay in a Gaussian way on the scale I. 



A. Ferromagnetic correlations at g = 

In contrast, the ferromagnetic spin-spin correlation 
function 



T n a n+RI 



«)(*■ 



n+Rl 



(«+r) (73) 



cannot be evaluated so easily. As is well known, in the 
ground state, C R Z can be written as a determinant of 
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an R x R Toeplitz matrix whose asymptote for large R 
can be obtained with the Szego limit theorem [3J|. Un- 
fortunately, in time-dependent problems the correlation 
function is not a determinant in general. However, below 
we avoid this problem in an interesting range of param- 
eters. 

Using the Jordan- Wigner transformation, Cjj z can be 
expressed as 

C Z R Z = {b axbia 2 ■ ■ .bR-ia R ) . (74) 
Here a n and b n are Majorana fermions defined 



(c* + c n ) and b r< 



c n . Using ([53)1 and (f54|) we get 



(75) 



{b m <l>n) = ' 
{(l>m<l>n) = ' 

(b m b n ) = ' 

The average in Eq. (|74|) is a determinant of a matrix when 
(a m a n ) = and (b m b n ) = for m ^ n, or equivalently 
when 3/3 m _„ = for m =/= n. Inspection of the last line 
in Eq. (|54|) shows that 9/3 m _„ « when |m — n| <C /. 
Consequently, when the correlation distance R <C Z we 
can neglect all 9/3 m _„ assuming that (a m a n ) = and 
{b m bn) — for m ^= n. In this regime, the correlation 
function is a determinant of the Toeplitz matrix 

[(Van+i>] m , n =i,...„R ■ (76) 

Asymptotic behavior of this Toeplitz determinant can be 
obtained using standard methods [34| with the result that 



exp 



-0.174 5 



cos 




(77) 



when 1 < i? < i. 

In this way we find that the final ferromagnetic corre- 
lation function at g = exhibits decaying oscillatory be- 
havior on length scales much less than the phase - ordered 
scale but both the wavelength of these oscillations and 
their exponentially decaying envelope are determined by 
£. As discussed in a similar situation by Cherng and Lev- 
itov [HI , this oscillatory behavior means that consecutive 
kinks are approximately anticorrelated - they keep more 
or less the same distance ~ £ from each other forming 
something similar to a ... -kink- antikink- kink- antikink-... 
lattice with a lattice constant ~ £. However, fluctuations 
in the length of bonds in this lattice are comparable to 
the average distance itself giving the exponential decay 
of the correlator Cjj z on the same scale of ~ £. 

We do not know the tail of the ferromagnetic correla- 
tion function when L 3> I because our approximations 
necessary to derive Eq. ([77]) do not work in this regime, 
but we can estimate that this tail is not negligible. In- 
deed, when R — I, then the envelope in Eq. (f77j) is 



exp 



I 

-0.174 - 



-0.174 In 



(78) 



Due to the smallness of the exponent 0.174 the tail is 
negligible only for extremely large tq. 



B. Correlations at the critical point 

In order to fill in the gaps in our analytic knowledge 
of the correlation functions, we attempted to make nu- 
merical simulations of the dynamical transition. As we 
wanted to get information on spin-spin correlation func- 
tions it proved convenient to work directly with spin de- 
grees of freedom rather than with the Jordan- Wigner 
fermions. We used the translationally invariant version of 
the real time Vidal algorithm [35[ . This algorithm, which 
is an elegant version of the density matrix renormaliza- 
tion group [36J, is an efficient way to simulate time evo- 
lution of an infinite translationally invariant spin chain. 
This ambitious task is made efficient by a clever trunca- 
tion of Schmidt decomposition between any two halves 
of the infinite spin chain. Our calculations of the en- 
tropy of entanglement demonstrate that for a finite tran- 
sition rate, the entropy saturates at a finite value beyond 
certain block size L. This saturation suggests that in 
our case the truncation of the Schmidt decomposition 
will make sense and, in principle, dynamical phase tran- 
sitions across quantum critical points can be efficiently 
simulated with the Vidal algorithm. 

In our simulations, we started the linear quench from 
the ground state at g = 10 which was prepared with the 
imaginary time version of the algorithm. The simulations 
were run for a range of tq such that the initial part of 
the evolution close to g = 10 was well in the adiabatic 
regime. We checked our results for convergence with re- 
spect to the truncation of the Schmidt decomposition (we 
used x U P to x = 40) and time step dt. We used fourth 
order Trotter decomposition. Wherever it was possible, 
we compared our numerical results with analytical re- 
sults which could be obtained for transversal magneti- 
zation, transversal spin-spin correlations, and ferromag- 
netic nearest-neighbor correlations. We also controlled if 
our truncation of the Schmidt decomposition is sufficient 
to preserve the norm of the state evolved in real time. 
As illustrated in panel A of Figure [31 our simulations 
were stable enough to cross the critical point and enter 
the ferromagnetic phase, but once in the ferromagnetic 
phase, the algorithm was breaking down. This is why we 
trust our numerical results at g — 1, but have no reliable 
results below g = 1. We can verify KZM at the critical 
point, but we cannot reliably follow the phase ordering 
in the ferromagnetic phase. 

In panel B of Figure[31 we plot the transverse spin-spin 
correlation Cf^ at g = 1 for several values of tq. For 
each tq, we plot both numerical correlator and its ana- 
lytic counterpart from Eq. (|70[) and they seem to approxi- 
mately coincide. Equation ([70)) can be also used to obtain 
analytically, but with some numerical integration, the ex- 
ponential tail of the transverse correlator when tq ^> 1: 



Cxx 
R 



0.44 / „ R 
exp -2.03— 



(79) 



accurate when R ^> £. This tail decays on the KZ cor- 
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FIG. 3: Panel A shows the dynamical transverse correlation 
Cf x as a function of magnetic field g in the linear quench. 
For each tq, we show both numerical (dashed) and analyt- 
ical (solid) result. The plots overlap near the critical point 
at g = 1 but diverge in the ferromagnetic phase when g < 1 
indicating a breakdown of our numerical simulations in this 
regime. Panel B shows analytic and numerical results for 
the dynamical transverse correlation function at the moment 
when the quench crosses the critical point at g — 1. The trans- 
verse correlators overlap well confirming that our numerical 
simulations are still accurate at the critical point. Finally, 
in panel C, we show the dynamical ferromagnetic correlation 
function CJj 2 at g = 1 and in panel D, we show the same cor- 
relation function after rescaling R/^/tq. The rescaled plots 
overlap quite well supporting the idea that near the critical 
point the KZ correlation length £ = ^/tq is the only relevant 
scale of length. 



relation length £ which proves to be the relevant scale of 
length. 

Encouraged by the agreement in transverse correla- 
tions in panel C we show the ferromagnetic spin-spin 
correlation functions at g = 1 for the same values of 
tq. They are roughly exponential and their correlation 

length seems to be set by £ = ^Jtq . To verify this scaling 
hypothesis we show in panel D the same plots as in panel 
C but with R rescaled as i?/£. We find the rescaled plots 
to overlap reasonably well confirming the expected s/tq 
scaling. The overlap is not perfect, but the scaling is ex- 
pected when tq > 1 which is not quite satisfied by the 
tq available from our numerical simulations. 



V. CONCLUSION 

Putting our analytical results and numerical evidence 
together, we are led to conclude that, in a quantum 



phase transition, the system initially follows adiabatically 
its instantaneous ground state. This adiabatic behavior 
becomes impossible sufficiently near the critical point: 
When crossing the critical regime the system gets ex- 
cited in a manner consistent with KZM, and imprinted 
with the characteristic KZ dynamical correlation length 
£ = Jtq. We find evidence for this correlation length 
both in correlation functions and in the entropy of en- 
tanglement - they are all determined by the same single 
length scale £. 

Once the system is excited, the non-trivial dispersion 
relation of its quasiparticle excitations leads to a grad- 
ual quantum phase ordering: Thanks to this post-critical 
evolution, the state of the system develops the second, 
longer, phase-ordered length scale which finally at g = 
becomes I = ^/tq In tq. This process makes short range 
ferromagnetic correlation function oscillatory rather than 
purely exponential, which means that on length scales 
shorter than I the random -kink-antikink-kink-antikink- 
train looks more like a regular crystal lattice. At the same 
time, thanks to phase ordering, a longer block of spins is 
necessary to saturate the entropy of entanglement. 

It is important to note that the first process depends 
on the universal characteristics of (quantum or classi- 
cal) second order phase transitions and can be modeled 
by KZM. Therefore, we expect that conclusions we have 
reached for the specific case of the quantum Ising model 
are generally applicable: Once the universality class of 
the transition is characterized by means of the relevant 
critical exponents, predictions of e.g. the entanglement 
entropy left in the wake of the phase transition can be 
made. By contrast, the dynamics of the phase ordering 
that follows can be model-specific, and is unlikely to be 
captured by the scalings of relaxation time and healing 
length that suffice for KZM. 
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